{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "5d904dee",
   "metadata": {},
   "source": [
    "# Example 6: Solving Partial Differential Equation (PDE)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d568912",
   "metadata": {},
   "source": [
    "We aim to solve a 2D poisson equation $\\nabla^2 f(x,y) = -2\\pi^2{\\rm sin}(\\pi x){\\rm sin}(\\pi y)$, with boundary condition $f(-1,y)=f(1,y)=f(x,-1)=f(x,1)=0$. The ground truth solution is $f(x,y)={\\rm sin}(\\pi x){\\rm sin}(\\pi y)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "0e2bc449",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "pde loss: 2.52e+00 | bc loss: 1.57e-03 | l2: 3.10e-03 : 100%|███████| 20/20 [00:25<00:00,  1.25s/it]\n"
     ]
    }
   ],
   "source": [
    "from kan import *\n",
    "import matplotlib.pyplot as plt\n",
    "from torch import autograd\n",
    "from tqdm import tqdm\n",
    "\n",
    "dim = 2\n",
    "np_i = 21 # number of interior points (along each dimension)\n",
    "np_b = 21 # number of boundary points (along each dimension)\n",
    "ranges = [-1, 1]\n",
    "\n",
    "model = KAN(width=[2,2,1], grid=5, k=3, seed=3)\n",
    "\n",
    "def batch_jacobian(func, x, create_graph=False):\n",
    "    # x in shape (Batch, Length)\n",
    "    def _func_sum(x):\n",
    "        return func(x).sum(dim=0)\n",
    "    return autograd.functional.jacobian(_func_sum, x, create_graph=create_graph).permute(1,0,2)\n",
    "\n",
    "# define solution\n",
    "sol_fun = lambda x: torch.sin(torch.pi*x[:,[0]])*torch.sin(torch.pi*x[:,[1]])\n",
    "source_fun = lambda x: -2*torch.pi**2 * torch.sin(torch.pi*x[:,[0]])*torch.sin(torch.pi*x[:,[1]])\n",
    "\n",
    "# interior\n",
    "sampling_mode = 'random' # 'radnom' or 'mesh'\n",
    "\n",
    "x_mesh = torch.linspace(ranges[0],ranges[1],steps=np_i)\n",
    "y_mesh = torch.linspace(ranges[0],ranges[1],steps=np_i)\n",
    "X, Y = torch.meshgrid(x_mesh, y_mesh, indexing=\"ij\")\n",
    "if sampling_mode == 'mesh':\n",
    "    #mesh\n",
    "    x_i = torch.stack([X.reshape(-1,), Y.reshape(-1,)]).permute(1,0)\n",
    "else:\n",
    "    #random\n",
    "    x_i = torch.rand((np_i**2,2))*2-1\n",
    "\n",
    "# boundary, 4 sides\n",
    "helper = lambda X, Y: torch.stack([X.reshape(-1,), Y.reshape(-1,)]).permute(1,0)\n",
    "xb1 = helper(X[0], Y[0])\n",
    "xb2 = helper(X[-1], Y[0])\n",
    "xb3 = helper(X[:,0], Y[:,0])\n",
    "xb4 = helper(X[:,0], Y[:,-1])\n",
    "x_b = torch.cat([xb1, xb2, xb3, xb4], dim=0)\n",
    "\n",
    "steps = 20\n",
    "alpha = 0.01\n",
    "log = 1\n",
    "\n",
    "def train():\n",
    "    optimizer = LBFGS(model.parameters(), lr=1, history_size=10, line_search_fn=\"strong_wolfe\", tolerance_grad=1e-32, tolerance_change=1e-32, tolerance_ys=1e-32)\n",
    "\n",
    "    pbar = tqdm(range(steps), desc='description', ncols=100)\n",
    "\n",
    "    for _ in pbar:\n",
    "        def closure():\n",
    "            global pde_loss, bc_loss\n",
    "            optimizer.zero_grad()\n",
    "            # interior loss\n",
    "            sol = sol_fun(x_i)\n",
    "            sol_D1_fun = lambda x: batch_jacobian(model, x, create_graph=True)[:,0,:]\n",
    "            sol_D1 = sol_D1_fun(x_i)\n",
    "            sol_D2 = batch_jacobian(sol_D1_fun, x_i, create_graph=True)[:,:,:]\n",
    "            lap = torch.sum(torch.diagonal(sol_D2, dim1=1, dim2=2), dim=1, keepdim=True)\n",
    "            source = source_fun(x_i)\n",
    "            pde_loss = torch.mean((lap - source)**2)\n",
    "\n",
    "            # boundary loss\n",
    "            bc_true = sol_fun(x_b)\n",
    "            bc_pred = model(x_b)\n",
    "            bc_loss = torch.mean((bc_pred-bc_true)**2)\n",
    "\n",
    "            loss = alpha * pde_loss + bc_loss\n",
    "            loss.backward()\n",
    "            return loss\n",
    "\n",
    "        if _ % 5 == 0 and _ < 50:\n",
    "            model.update_grid_from_samples(x_i)\n",
    "\n",
    "        optimizer.step(closure)\n",
    "        sol = sol_fun(x_i)\n",
    "        loss = alpha * pde_loss + bc_loss\n",
    "        l2 = torch.mean((model(x_i) - sol)**2)\n",
    "\n",
    "        if _ % log == 0:\n",
    "            pbar.set_description(\"pde loss: %.2e | bc loss: %.2e | l2: %.2e \" % (pde_loss.cpu().detach().numpy(), bc_loss.cpu().detach().numpy(), l2.detach().numpy()))\n",
    "\n",
    "train()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e2246bab",
   "metadata": {},
   "source": [
    "Plot the trained KAN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "02e2a0ba",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAFICAYAAACcDrP3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABRrUlEQVR4nO3dd1gU1/oH8O8sS+8goggWEJWiBBW7IolKNOWqydXYNXaNGk2uvWAJ9tg1iUajaCS5aoqxRb1gr1hQQEURFRDpZVl2YXfO7w+z8wMrZdnZhffzPDxP4rI7Lzt75rtzzpkzHGOMgRBCCNEiidgFEEIIqX4oXAghhGgdhQshhBCto3AhhBCidRQuhBBCtI7ChRBCiNZRuBBCCNE6ChdCCCFaR+FCCCFE6yhcCCGEaB2FCyGEEK2jcCGEEKJ1FC6EEEK0jsKFEEKI1lG4EEII0Tqp2AUQYggYY8jMzIRMJoOVlRUcHR3BcZzYZRGit+jMhZA3yMnJwbp16+Dp6QknJyc0atQITk5O8PT0xLp165CTkyN2iYToJY7uREnIqx07dgyffPIJ5HI5gOdnLxqasxYLCwvs378fwcHBotRIiL6icCHkFY4dO4YPPvgAjDHwPP/a35NIJOA4DocOHaKAIaQEChdCXpCTkwNXV1cUFha+MVg0JBIJzM3NkZSUBDs7u6ovkBADQGMuhLxg586dkMvlZQoWAOB5HnK5HLt27ariyggxHHTmQkgJjDF4enoiISEB5WkaHMfB3d0d8fHxNIuMEFC4EFJKRkYGnJycKvV8R0dHLVZEiGGibjFCSpDJZJV6fn5+vpYqIcSwUbgQUoKVlVWlnm9tba2lSggxbBQuhJTg6OgIDw+Pco+bcBwHDw8PODg4VFFlhBgWChdCSuA4DpMmTarQcydPnkyD+YT8gwb0CXkBXedCSOXRmQshL7Czs8P+/fvBcRwkkjc3Ec0V+gcOHKBgIaQEChdCXiE4OBiHDh2Cubk5OI57qbtL82/m5uY4fPgwevToIVKlhOgnChdCXiM4OBhJSUlYu3Yt3N3dSz3m7u6OtWvXIjk5mYKFkFegMRdCyoAxhoiICLz33ns4efIkgoKCaPCekDegMxdCyoDjOGFMxc7OjoKFkLegcCGEEKJ1FC6EEEK0jsKFEEKI1lG4EEII0ToKF0IIIVpH4UIIIUTrKFwIIYRoHYULIYQQraNwIYQQonUULoQQQrSOwoUQQojWUbgQQgjROgoXQgghWkfhQgghROsoXAghhGgdhQshhBCto3Ah5C2Ki4uRnJyMuLg4AMCDBw+QlZUFnudFrowQ/UW3OSbkNXJycrB//37s2bMHMTExyM/PR1FREczMzODk5ITOnTtj5MiR6NixI6RSqdjlEqJXKFwIeYULFy5g6tSpiI6ORkBAAD744AO0aNECVlZWyMnJQVRUFA4ePIj79++jf//+WLJkCZycnMQumxC9QeFCyAv+/vtvDB8+HFZWVli6dCl69eqFoqIihIeHQ6lUwsbGBp999hmKi4sRHh6OkJAQ+Pj4ICwsDM7OzmKXT4heoHAhpIR79+7h/fffh6WlJcLDw+Ht7Q2O45CQkICWLVsiNzcXjRo1QlRUFOzt7cEYw9mzZzFw4EB07doV27Ztg6mpqdh/BiGiowF9Qv6hVqsRGhqK7OxsbNy4UQiWN+E4Dp06dcKKFSvwxx9/4OjRozqqlhD9RuFCyD/u37+PgwcPom/fvujUqdNbg0WD4zj07t0b7dq1w9atW6FSqaq4UkL0H01xIeQf58+fh0wmwyeffILExEQUFBQIjyUlJUGtVgMAioqKEBMTAxsbG+FxFxcX9O3bFyEhIUhNTYWrq6vO6ydEn1C4EPKPO3fuwMLCAu7u7hg7dizOnTsnPMYYg1KpBACkpKSge/fuwmMcx2H16tVo3rw55HI5UlJSKFxIjUfhQsg/CgsLIZVKYWpqCqVSCYVC8crfY4y99JhKpYK5uXmpECKkJqNwITXew4cPERERgdOnT0MulyMnJwdt27aFpaWl8DuFhYU4f/68ECIdOnQQLpzkOA7169dHWloaVCoV4uPjERAQADMzM7H+JEJER1ORSY3z+PFjREZGIiIiAhEREXj06BE4jkOjRo3w6NEjbNq0CaNGjSr1nISEBAQEBCA3NxcNGzbE1atXYWdnJzzOcRxmz56NlStXgud5mJqaon379ujatSuCgoLQtm1bmqJMahQKF1LtJScnIyIiQgiUhIQEAICfn59w8O/SpQt4nkenTp1gb2+Po0ePlhqwf911LsDzbrKUlBQEBgbio48+wrBhwxAZGYnIyEicOnUKOTk5wtmOZnsBAQEwMTER5f0gRBcoXEi1k5qaWipM4uPjAQC+vr7CwT0wMBCOjo4vPXfTpk346quvMHfuXMycOVPo+npTuCgUCnz55Zc4ePAg/ve//6Fp06bC66nVaty8eVOo5fTp08jLy4OFhQU6duyIoKAgBAUFoVWrVjA2NtbBu0OIblC4EIOXlpYmnClERETgzp07AAAvL69SYVK7du23vlZBQQE+//xzHD58GAsXLsT48eNhZmaGhw8fok2bNkK32OXLl2FnZ4f8/Hx88803+P7777FmzRqMGDHija+vUqlw/fp1IfzOnDkDmUwGKysrdOrUSQgbf39/WgyTGDQKF2JwMjIycOrUKSFMYmJiAABNmjQRwqRr166oU6dOhV4/PT0dEydOxF9//YXg4GBMnToVXl5euHv3Lnieh4mJCRo3bozLly9j1apVuHHjBhYtWoTx48fDyMioXNsqLi5GVFSUEDZnz56FXC6HjY0NOnfuLISNn59fuV+bEDFRuBC9l52djVOnTgkH4OjoaACAh4dHqTCpV6+e1rZZUFCArVu3Yv369Xj27Bnc3d3h6ekJa2trZGdn4+7du0hJSUGrVq2wYMECBAYGQiKp/IIXRUVFuHLlihCc586dg0KhgJ2dHbp06SKETfPmzbWyPUKqCoUL0Tu5ubk4ffq0ECY3btwAYwwNGzYUgiQoKAhubm5VXktqaipOnjyJU6dOISEhAQqFAvb29vD19UWPHj3Qtm1bWFhYVNn2lUolLl26JITNhQsXoFQq4eDggMDAQCFsfHx8yrxcDSG6QOFCRJefn48zZ84IB9Br166B53m4uroKB8+goCA0bNhQ1DrVajUYY5BIJKKdNSgUCly4cEF4ry5evIji4mI4OTmVCptmzZpR2BBRUbgQnZPJZDh37pxwZnL16lWo1Wq4uLiUOjNxd3enA+RbyOVynD9/Xgiby5cvQ6VSwdnZWXgfg4KC4OnpSe8l0SkKF1LlNAdATZiUPACWDBM6AFaeTCYT3uuIiIhSwV0ybCi4SVWjcCFap+m60YRJya6brl27Cgc56rqpenl5ecJZYskuRzc3t1JhI3aXI6l+KFxIpWkGnTVhUnLQuWSY0KCz+HJyckqNb5WcLFEybHQxWYJUbxQupNw002U134bPnz8vTJcNDAwUDlI0XVb/ZWVl4cyZM8K+1Ezzdnd3LzWZwsXFReRKiaGhcCFvVfJCP821F5oL/TTXXnTt2pUu9KsGMjIyhGngJS9Q9fT0FIKmMheokpqDwoW8pOQSJRERETh79qywREnnzp2FMxNaoqT6S0tLEy5gLbm0TrNmzUqFjZOTk8iVEn1D4UKExRU1B5AzZ84Iiytq1rvq2rUrLa5I8PTp01Jho1kU1MfHRwib1y0KSmoWCpcaiOd53Lp1SzhAnD59Gjk5OTAzMxNW6u3atSstC0/eKjk5udS9cTS3M2jRooUQNl26dBFWkCY1B4VLDcAYQ0xMjHAAOHXqFLKysoQbWmnChG5oRSrrdTdie+edd4Sw6dy5M2xtbcUulVQxCpdqiDGGO3fulAqT9PR0GBsbo127dkKYtG/fnm7FS6rUw4cPS4VNUlISJBIJWrZsKYRNp06dYG1tLXapRMsoXKoBxhji4+OFBhwZGYlnz55BKpWiTZs2Qph06NChShdZJORNGGNISEgQPqcRERF4+vQpjIyM0Lp1ayFsOnbsCEtLS7HLJZVE4WKAXmykkZGRSElJeamRdujQAVZWVmKXS8grleVLUVBQENq3b09figwQhYuBSExMLNUInzx58lL3QseOHUvd950QQ/Jid25kZCQyMjJgYmKCtm3bUneugaFwMRDNmjXDvXv3Sg2MdurUCXZ2dmKXRkiV4HkesbGxL01ECQ8PR//+/cUuj7wFhYuB0OwmWpuL1FQlD1XUDvQfhQshhBCto7U7tEQzOJmZmSl2KZUikUjg6+tLs3VIuVEbICXRmYuW8DyPiRMnwtXVVWsztNRqNTiO0+nKwmfOnMH8+fPRokULnW2TVA9V0QaqmkqlgpGRUaluNmoD2kFnLlpkamqKUaNGwdnZuVKvw/M8jh07hk2bNsHGxgYzZ85E8+bNq6yfWaFQIDs7G87OzpDJZKDvG6SiytsGGGPIycnB/fv34ezsDFdXV518mSouLsauXbsQHh4OLy8vzJ49W6iZ2oB2ULjoGcYYIiMjMWTIEGRlZQEAYmJicOjQIbi6ulbJ9o4fP45Jkyahb9++1BVAdIYxhlu3bmHUqFG4desWHBwcEBISgs8//7xKb93AGENYWBgmTZoEhUKBkydPIi0tDT/99BMtf6RFdCcnPSOTybBgwQJkZWXB398fLi4uuHXrFrZs2QKe57W+PZVKhR07duDx48dITk6mJfSJzuTm5mLSpEm4evUqJBIJUlNT8fXXX+Pvv/+u0jOH5ORkhIaGQqFQoF27djA3N8cff/yBs2fPVtk2ayIKFz2iOYu4dOkSbG1tsWXLFkybNg0cx+Hnn39Gamqq1repuWjNxMQEw4YNo3AhOsEYw759+3D+/Hk4ODggPDwcPXv2RH5+PkJCQpCbm1tl2w0PD8fDhw/RsGFD7Nq1C127doVSqcS+ffuoO0yLKFz0iEqlQlhYGFQqFXr16oWWLVuif//+cHV1xZMnT7T+jY4xhmPHjiE3NxdeXl7o1KmT1l6bkDeRy+X48ccfoVarMXDgQPTq1QtLliyBo6Mjrl27hj///LNKDvQymQzh4eFgjGHw4MHw8PBA3759wXEcTpw4YfAz3fQJhYseSUhIwJkzZ2BsbIxBgwbByMgIdevWRY8ePcDzPP744w+o1Wqtba+oqAhHjhwBAAQHB9PKtEQnGGO4du0abt68CWtrawwbNgwSiQTNmzfHJ598ArVajR07dkChUGh9u9evX0dMTAxsbGzw6aefguM4vPvuu3BycsLjx48RFRWl1W3WZBQueoIxhpMnTyI7OxseHh5o3749OI4Dx3H48MMPYWRkhEuXLuHp06da22ZSUhKio6NhYmKC4OBgrb0uIW/z559/QqFQoFWrVvDx8RGm3A8bNgwWFha4cuUKbty4ofWzl0OHDkGpVKJly5Zo2rQpAKBevXpo3rw5VCoVLly4oNXt1WQULnpCpVLh8OHDYIyhW7duwpphHMehbdu2cHFxQVpaGi5duqSVBscYQ1RUFLKzs+Hi4gJfX99KvyYhZSGTyXD8+HEAwEcffSTM0NLcVKx169aQy+X47bfftLrdgoICnDx5EgDQs2dP4S6rxsbGaN++PQDg4sWLUKlUWt1uTUXhoieePn2KqKgoSKVSvP/++6Uec3JyQuvWrcHzPCIjI7W2zcjISPA8j5YtW8LBwUFrr0vIm9y5cwf379+HlZUV3nvvvVLXb5mZmaF3794AgKNHjyIvL09r233w4AHu3bsHCwsLvPvuu8J2OY5D9+7dMWLECAwdOlRr26vpKFz0gOYsIiMjA3Xq1EHLli1LNTiJRIKgoCAAwPnz51FQUFDpbRYWFuLKlSsAgC5duuh0FQBSczHGcObMGcjlcjRt2hQeHh6lHuc4Dj169ICtrS3i4+MRHR2tte2ePXsWMpkMHh4eaNKkSanHO3XqhG3btmHgwIE0Y1JL6IiiJyIjI6FWq+Hv749atWqVeozjOOGGSQ8ePMCjR48qvb2UlBQkJCTA1NQUrVu3plVmiU6o1WqcOnUKwPMD+qtuAubu7o4WLVpAqVTi+PHjWukG5nle2O6rbqKnGd8k2kPhogcKCwuFgcTAwMBXXp3s4eEBNzc3yGQyREVFVbrB3b59G7m5uXB2dkbjxo0r9VqElFVmZiZu3LgBIyMjBAYGvvJ3TE1N0a1bNwDAyZMntTJrLCcnB9euXYNEIkGXLl0q/Xrk7Shc9EBSUhLu378PMzMzdOjQ4ZXfoGxsbODv7w/GGM6fP1+p7THGcPnyZfA8D29vbxpvIToTGxuL1NRU2Nvb45133nnlZ10zPdjU1BSxsbF4/PhxpbcbHx+P5ORk2NjYoFWrVnSWogMULiLTzL3Py8uDi4sLPD09X/l7mq4xAIiKikJhYWGFt6lSqXD58mUAQJs2bap0HSdCNBhjuHjxIoqLi9G0aVPUrVv3tb/r7e2NBg0aIC8vDxcvXqzUmbrmy5RSqYSHh0eVrNFHXkbhogcuXLgAnufh5+f32tsWcxyH1q1bw9TUFAkJCZW63iUzMxN3796FkZER2rRpQ9/iiE6o1WpcvHgRwPMvNZqpwK9ia2uLNm3aCAu5ViZceJ4Xup0DAgJgbm5e4dciZUfhIjKlUinM2mrfvv0bZ215eHjA2dkZeXl5uH37doW3mZCQgPT0dNjZ2cHLy6vCr0NIeeTk5OD27duQSCTCWfjrcByHwMBAcByHy5cvIz8/v8LblclkuHnzZqmzf1L1KFxElpqaivj4eJiYmCAgIOCNZxH29vZo1qwZ1Go1rly5UqFvc4wx3Lx5E0VFRWjQoEGl7z1DSFlpzritra3fen8ijuPQpk0bWFlZ4dGjR0hISKjwdh89eoQnT57A0tLyteM8RPsoXEQWFxeH7OxsODk5CctRvI5UKkVAQAAA4OrVqxW+klizflLz5s1hZmZWodcgpDw064kpFAo0bNgQ9erVe+tzGjZsiIYNG0Iul+Py5csV/jJ1/fp1yOVyuLm5oX79+hUpn1QAhYuIGGO4cuUK1Go1mjZt+tZZW5pxF4lEgjt37iAnJ6fc21QoFMKFaTRrhuiKZlAdAN55551XXt/yIktLS+HL1NmzZys87qIJphYtWtDirDpE4SIiTfcWALRu3RrGxsZvfY63tzdsbGyQlpZWoa6CtLQ0JCYmwsTEBH5+fhQuRCcUCgVu3rwJAG/t/tXgOA4dO3YEx3G4fv06ZDJZuberVCpx7dq1cm2XaAeFi4hyc3MRGxsLiUQifEN7GxcXFzRo0AAKhaJCq8beu3cPOTk5cHR0fGnpDUKqytOnT5GYmAhTU1P4+/uXOVxatmwJCwsLPH78GImJieXe7rNnz5CQkABjY2M6U9cxChcRPXz4EKmpqbC2toavr2+ZPvjm5ubw8/MDAOGsp6w0/c8qlQoeHh5wdHSsUN2ElAdjDLGxscjLy4Ozs3O5vtQ0bNgQ9evXR0FBAa5du1buL1OaMc1atWq9tJ4YqVoULiLRHOgLCwtRv379Mg1wAv8/7gIAN27cKNfFlJpBVeB5v/ebrjMgRJuuXr0qjC3a29uX+XnW1tbw8/MTLsAsD8aYMPGlSZMm9GVKxyhcRKQZ4PTz8yvTACfw/10FpqamSExMRGpqapm3V1BQgJiYGOE1CNEFtVqNq1evAng+tlieVYdLXpty9erVcn2Z4nm+3GOaRHsoXERSWFiIGzduACj/QKOHhwecnJyQm5uLuLi4Mj8vOTkZT548gZmZGVq0aEH9z0QncnJyEBcXB4lEUu4VuEuuTPHw4UMkJyeX+bkVGdMk2kPhIpLU1FQ8fPgQpqamL92/5W0cHBzQtGnTcl1MyRhDTEwMZDIZ6tSpgwYNGlSmfELKTDO2aGNjAx8fn3I/v3Hjxqhbty5yc3MRHR1d5nGXR48e4enTp7CysnrrRZtE+yhcRBITE4Pc3Fw4OTmVe8n7khdTaq6TKYuoqCjwPI9mzZrB1ta23DUTUl6asUWFQoEGDRrAxcWl3K9hb28PHx8f8Dxf5nEXxpgwJunm5kaLVYqAwkUEmoFGtVqNZs2aVWjJ+zZt2kAikSA2NhbZ2dlv/X2VSlWq35tWQia6wBjDpUuXAJRvbLEkIyMjtGvXDsDzccqioqIyPe/SpUtgjMHPzw+Wlpbl3i6pHAoXEZRc8j4gIKDct1XlOA6+vr6wtbVFamoq7t+//9bnZGVl4c6dOzAyMqI7TxKdKTm2WNEVuDmOQ9u2bSGVSnHv3j2kp6eXabuamZG08rc4KFxEkJ2djbi4OGHJ+4pwcXGBh4cHlEplmdZdunfvHtLS0mBnZ1ehfm9CKiIlJQUPHz6EmZlZpS5i9PLyQq1atZCZmYnY2Ni3/v7Tp0/pNt4io3ARwYMHD/Ds2TPY2NiU+eLJF5mbmwvjLufPnwfP86/9XcYYLly4gKKiImFwlJCqxhhDdHQ08vLyUKdOnUqtCOHk5AQvLy+oVCpcuHDhjV+mGGN0G289QOGiY5rxFqVSCXd390od6Euuu5Sbm/va31Or1Th16pTwHFoJmejKxYsXwfM8fH19X3sjvLIwNjYWxl0uXLjw1hXBL168CLVaDR8fH7qNt0goXHSMMYbz588DeD6wXtG74nEch1atWsHW1hZJSUm4c+fOa383LS0NN2/ehJGREQIDAyu0PULKq6ioSJjd1a5du0pNIuE4Dh06dICRkRFu376NjIwMnWyXVByFi47l5eXhxo0bQmOpjPr166Np06ZQKBQ4ffr0a7sKoqOjkZaWBicnJ7pZEtGZlJQU3LlzByYmJujQoUOlP3fNmzeHo6Mj0tPT33gn1mfPniEuLg7GxsZo3749fd5FQuGiZWq1GtnZ2UhMTHzlOMiDBw/w5MkTWFlZlfviyReZmZkJZyInT55EcXHxS7/DGMOpU6egUqnQokUL1KlTp8LbI6SsNNeZZGdno06dOlq5nXadOnXg4+OD4uLi197fhTGGW7duITMzE05OTvD29q70dknFULho2YULF9CmTRt89NFHL93MS7P4nlwuh4eHh1auku/evTtMTExw/fp1PHr06KXHCwsLERERAQAIDAws97RnQirq9OnT4Hkefn5+Wlk00tjYGJ07dxZe+3XXu5w+fRpqtRrNmzeHk5NTpbdLKobCRcscHByEG3IlJSWVeozneWFgvW3btpW+sIvjOPj7+6Nhw4bIzs7G8ePHS32b0yz5EhMTAwsLC3Tr1o26CIhOyOVynDlzBsDzLzXaGPfgOA6BgYEwNjbG7du3X7nOmEKhELbbpUsX+jIlIgoXLXNxcUHdunUhl8sRGxtb6mCfmZmJK1euQCKRICgoSCvbs7Ozw/vvvw/GGPbt2weFQlHq8f3796OgoAAtWrSg61uIzty9exd3796FhYUFAgMDtfalxtfXF/Xq1UN2djYuXrz4UtdYYmIi4uLihC5j+jIlHgoXLbO2tkbTpk2FPmcNzf+npKTAwcFBa1cNcxyHf//737CwsMDly5dLXVD57Nkz7N+/HxzHoV+/fhWemUZIeTDG8Pfff6OgoABeXl5o2rSp1l7bwcEB7du3B8/zOHr0aKlxTcYYTp8+jby8PDRs2JDGW0RG4aJlEokE/v7+AJ4vFFlykP3YsWMoLi6Gv79/hRbwexXNlOROnTpBLpdjw4YNUCqVYIwhLCwMCQkJqFevHvr27Uvf4ohOKJVKHDlyBAAQHBys1XW9JBIJevbsCYlEgjNnzpRaCkalUuHw4cNgjCEwMBA2NjZa2y4pPwoXLeM4Du3atYNUKkVsbKzw4c/NzcWxY8cAAL169dLqjYtMTU0xefJkmJmZ4dChQ9ixYwciIyOxZs0aAMCIESNoVViiM/fv38fNmzdhZmaGnj17avVLDcdx6NKlC5ydnZGUlFRqCv6jR49w4cIFSKVSfPjhh/RlSmQULlWgefPmqF27NtLS0nD16lXhwsn4+HjY29ujR48eWm9w3bp1w6BBg1BUVIRp06bhX//6F1JTU9GqVStMmDABEgntaqIb2dnZ8PDwQPPmzdG8eXOtv369evXQtWtXqNVqhIeHQ6VSgTGGgwcPIjMzE+7u7mjXrh2Fi8joiFMFnJ2d0bZtW6jVavz2229QKpXYsWMHiouL0aVLlypZ68jY2BihoaHo378/JBIJVCoV2rZti61bt6J27dpa3x4hr9OpUyecOnUK+/btg7W1tdZfXyKRYODAgTAxMcH//vc/3Lx5Ezk5OQgLCwNjDH379qUlX/QAzdOrAkZGRujXrx/+/PNPHDp0CBs3bsSRI0dgamqKkSNHVsn0SI7jUKtWLWzfvh23b9+GRCKBh4cHbGxs6Bsc0SmO42BlZQUrK6sqe/3AwEAEBATg3LlzCA0Nhb+/P27duoVatWph8ODB9JnXAxQuVYDjOPTo0QPvvPMOoqKiMHv2bKhUKnz00Ud47733quyDz3EczMzM0Lp16yp5fUL0haWlJWbMmIHPPvsMf/75J/766y/wPI/Ro0drdXYaqTjqFqsidnZ2WLVqFTw8PCCRSNCuXTusXLmSViQmRAs4jkNwcDDmzZsHGxsbSKVS9OvXD19//TUtVKkn6MxFixhjyM7OFmaC+fr64rfffsOTJ0/QpEkT2Nvbl+mWxGJSKpVil0AM2IttoKp9/vnn6Ny5M+RyOby8vMAYQ1ZWVqVek9qAdlC4aAnHcahfvz42bNjwym9Omnn/+q6wsBC2trZil0EM0NvaQFX7888/tfI61Aa0g2Nvuz8uKRPG2FtvNWwoOI6jAVFSbtQGSEkULoQQQrSOBvQJIYRoHY25GIiSJ5h0uk5qKmoHhoPOXAzE9evXIZFIcP36dbFLIUQ01A4MB4ULIYQQraNwIYQQonUULoQQQrSOwoUQQojWUbgQQgjROgoXQgghWkfhQgghROsoXAghhGgdhQshhBCto3AhhBCidRQuhBBCtI7ChRBCiNZRuBBCCNE6ChdCCCFaR+FCCCFE6yhcDABjDNnZ2QCA7OzsanOfckLKg9qBYaFw0WM5OTlYt24dPD090a1bNwBAt27d4OnpiXXr1iEnJ0fcAgnRAWoHholjFP966dixY/jkk08gl8sBvPr2rhYWFti/fz+Cg4NFqZGQqkbtwHBRuOihY8eO4YMPPgBjDDzPv/b3JBIJOI7DoUOHqGGRaofagWGjcNEzOTk5cHV1RWFh4RsblIZEIoG5uTmSkpJgZ2dX9QUSogPUDgwfjbnomZ07d0Iul5epQQEAz/OQy+XYtWtXFVdGiO5QOzB8dOaiRxhj8PT0REJCQrlmwnAcB3d3d8THxwv90IQYKmoH1QOFix7JyMiAk5NTpZ7v6OioxYoI0T1qB9UDdYvpEZlMVqnn5+fna6kSQsRD7aB6oHDRI1ZWVpV6vrW1tZYqIUQ81A6qBwoXPeLo6AgPD49y9xdzHAcPDw84ODhUUWWE6A61g+qBwkWPcByHSZMmVei5kydPpkFMUi1QO6geaEBfz9D8fkKoHVQHdOaiZ+zs7LB//35wHAeJ5M27R3Nl8oEDB6hBkWqF2oHho3DRQ8HBwTh06BDMzc3BcdxLp/mafzM3N8fhw4fRo0cPkSolpOpQOzBsFC56Kjg4GElJSVi7di3c3d1LPebu7o61a9ciOTmZGhSp1qgdGC4aczEAjDFERETgvffew8mTJxEUFESDlqTGoXZgWOjMxQBwHCf0JdvZ2VGDIjUStQPDQuFCCCFE6yhcCCGEaB2FCyGEEK2jcCGEEKJ1FC6EEEK0jsKFEEKI1lG4EEII0ToKF0IIIVpH4UIIIUTrKFwIIYRoHYULIYQQraNwIYQQonUULoQQQrSOwoUQQojWUbgQQgjROgoXQgghWkfhoudkMhnu3buHW7duAQBSU1NRVFQkclWE6FZxcTGSk5MRFxcHAHjw4AGysrLA87zIlZHXodsc66mEhARs27YNf/75J548eYLi4mIolUrY2NjA398fw4YNQ9++fWFtbS12qYRUmZycHOzfvx979uxBTEwM8vPzUVRUBDMzMzg5OaFz584YOXIkOnbsCKlUKna5pAQKFz2jVquxd+9ezJ49G4WFhejZsye6d++O+vXrg+d53L9/H0eOHEFERARatmyJDRs2wNvbW+yyCdG6CxcuYOrUqYiOjkZAQAA++OADtGjRAlZWVsjJyUFUVBQOHjyI+/fvo3///liyZAmcnJzELpv8g8JFj/A8j++++w7Tp09Hly5dsGzZMvj4+ODy5cu4cuUKAKBHjx7w8PDAhQsX8NVXXyE/Px/79u2Dr6+vyNUToj1///03hg8fDisrKyxduhS9evVCUVERwsPDhTP4zz77DMXFxQgPD0dISAh8fHwQFhYGZ2dnscsnAMCI3oiIiGB2dnbs008/ZVlZWYznecYYY3PnzmUAGAAWFhbGGGOM53n26NEj1qFDB9apUyeWnZ0tYuWEaM/du3dZo0aNmK+vL7t9+7bQDh48eMBsbW0ZANaoUSOWlZXFGHveFk6fPs1cXV3Z4MGDmUKhELN88g8a0NcThYWFWLRoEZydnbFmzRrY2dmB47jX/j7HcXBzc8OGDRtw79497N69W4fVElI11Go1QkNDkZ2djY0bN8Lb2/uN7QB43hY6deqEFStW4I8//sDRo0d1VC15EwoXPREVFYWLFy9iwoQJqFev3lsbFPC8Ub3zzjvo168ffvrpJ8jlch1USkjVuX//Pg4ePIi+ffuiU6dOZWoHwPO20Lt3b7Rr1w5bt26FSqWq4krJ29D0Cj0RGRkJU1NTdOvWDXFxcaUax7Nnz4T/fvz4MaKjo4X/t7OzQ+/evbF7924kJibS4D4xaOfPn4dMJsMnn3yCxMREFBQUCI8lJSVBrVYDAIqKihATEwMbGxvhcRcXF/Tt2xchISFITU2Fq6urzusn/4/CRU/cvXsXtWvXhrGxMbp164a0tDThsZJBExISgsWLFwv/P2DAAMyfPx9SqRSPHz+mcCEG7c6dO7CwsIC7uzvGjh2Lc+fOCY8xxqBUKgEAKSkp6N69u/AYx3FYvXo1mjdvDrlcjpSUFAoXkVG46AHGGBQKBUxNTWFkZASFQgGFQvHK3y0uLkZxcbHw/0VFRTAxMYFEIsGVK1fQsmVL1K5dW1elE6IVRUVFuHfvHqKjoyGVSmFqagqlUvnadqBpMyWpVCqYm5uXCiEiHgoXPcBxHGrVqoXLly9DrVYjKCgIOTk5wuPx8fFISEgAADRv3hwuLi7CYy1atEBOTg5kMhnmz5+P+fPno1atWvD19YWPj0+pH0dHR13/aYSUolKpEB8fj5iYmFI/9+7dE87QLSwskJOTg7Zt28LS0lJ4bmFhIc6fPy+ESIcOHYQLJzmOQ/369ZGWlgaJRAJ7e3tR/j7y/yhc9ESrVq2wc+dOpKamYs+ePaUemz9/PkJDQwEAX3/9NQYNGiQ8xnEcwsLCYGNjg0OHDiEjI0NosBEREfj++++FRuvs7AwfH5+XgsfOzk5nfyepGdRqNR48ePBSiNy9e1dYvqhWrVrw8fFB165d8cUXX8DHxwdPnz7F0KFDcfnyZSxfvrzUayYkJCAgIAC5ublwdnbGL7/8Uuqzy3EcZs+ejTp16lCXmB6gcNEDaWlpuHjxIgoLC7Fz5060b9++1FIWEomk1H8bGRkJ/y+Xy7Fr1y506tQJnTt3hpGRET799FPhcU13Q8kGfuzYMWzcuFFYl8nFxeWl0PH29i41WErIq/A8j4cPH74UInFxcULXlL29PXx8fNChQweMGTNG+Iy9qvs2MzMT7u7u2LlzJwYMGFDqM1jyc89xXKm2wBhDSkoK/vvf/6J58+b02dUDFC4iysjIwMqVK7Fx40ZIJBK0bdsWv/76K/r06YNevXq9dRomz/P46aefcP36dfz++++lGp+GiYkJfH19X7qCX6FQ4O7du6UOCAcPHsTatWvB/lm0wc3N7ZWhU7KrgtQMPM/j8ePHrwwRzRR4Gxsb+Pj4ICAgAMOHDxc+M3Xq1CnzlGJHR0d88cUX+Oqrr7B+/XrMnDmzTGuGKZVKLFq0CElJSUhISED79u2xcOFCBAcHl3nbRLsoXESQmZmJ1atXY8OGDQCAqVOnYtq0aSgqKsLHH3+M8ePHY8eOHQgKCoJEIoFEIoFUKgXHceA4DowxqNVqhIeHY8GCBRg/fjw6duxYrhrMzMzg5+cHPz+/Uv8ul8tx586dUgeQ/fv3Y9WqVcLvNGzY8KXQ8fLygrm5eeXfHCIqxhiSkpJeCpHY2FjIZDIAgJWVFby9veHn54eBAwcKn4GyXp/1NsOHD8fp06exfPlyWFhYYPz48TAzMwMASKVSSKXSUmcs+fn5+Oabb/DLL79gy5YtqF+/PubPn4+ePXsKIdOtWzcKGR2jtcV0KDs7G99++y3WrVsHnucxadIkfPXVV6hVq5bwO7GxsRgyZAgSExMxfvx4jBgxAjzPIyUlBQDQqFEj5ObmYvPmzdi7dy8GDx6MFStWwMLCokprl8lkiIuLe+mg8/jxYwDPuync3d1fCp2mTZsKBwaiPxhjePr06Uv7MyYmBnl5eQCeD6x7eXm9NEZXv379Kj9Qp6enY+LEifjrr78QHByMqVOnwsvLC3fv3gXP8zAxMUHjxo1x+fJlrFq1Cjdu3MCiRYswfvx4GBkZgTGGv//+G/Pnz8fly5fRqVMnLFq0CEFBQVVaN/l/FC46kJOTg7Vr12LNmjUoLi7GF198gf/85z+vXcE1OTkZixcvxi+//AKpVApvb2+4ublBrVYjMTERd+/ehaOjI2bMmIEhQ4bA1NRUx3/R/8vLy0NsbOxLB6jk5GQAz8eIGjduLByYNAeqJk2awMTERLS6awrGGNLS0l4ZItnZ2QCen8U2a9bspRBp2LBhqfE+XSsoKMDWrVuxfv16PHv2DO7u7vD09IS1tTWys7Nx9+5dpKSkoFWrVliwYAECAwNfqpcxhiNHjmD+/PmIiopCYGAgFi5ciMDAQJH+qpqDwqUK5eXlYd26dfj222+hUCgwYcIETJ8+vUyrtqrVasTFxeHQoUO4fPky0tLSYGxsjEaNGiEoKAg9evTQ6+tZcnJyXnlAS01NBfC8e8PT0/Ol0GncuDGMjY1Frt4wlZwpqPm5ffs2MjMzATwff2vatOlLIeLu7v7K8Tp9kZqaipMnT+LUqVNISEiAQqGAvb09fH190aNHD7Rt2/atZ+6MMfz1119YsGABrl+/jnfffRcLFy5Ep06ddPRX1DwULlUgPz8fGzZswKpVqyCXyzFu3DjMmDEDdevWrdDracZYOI7T64NAWWRmZr4ydNLT0wEAxsbGaNq06Uuh4+HhYfB/u7ZkZ2e/FCAxMTHCqg5SqbTUe6j5ady4scHfUEutVoMxJoxFlhdjDH/88QcWLFiA6OhodO/eHQsXLkT79u2roNqajcJFi2QyGTZt2oSVK1ciPz8fY8aMwcyZM1GvXj2xS9N7r+u6ycrKAgCYmpqiWbNmL4VOo0aNRO26qUq5ubmluhw1IfL06VMAz6fmljz70/x4enpSl+Nb8DyP3377DQsWLEBMTAzef/99LFy4EG3atBG7tGqDwkUL5HI5Nm/ejOXLlyM3NxcjR47E7Nmz4ebmJnZpBo0xhmfPngkH1ZI/ubm5AABzc3N4eXm9FDr169c3mNCRyWSvDJGkpCQAzydLlBy3KjlZQszxtuqA53ns27cPISEhiIuLwwcffICFCxeiVatWYpdm8ChcKqGwsBDfffcdli1bhqysLIwYMQJz5sxBgwYNxC6tWtNcMPeq0NFMl7W0tIS3t3epg7Gvry9cXV1Fm5Iql8tLzbjT1P/o0SPhdzQz7kr+NGvWjKZ5VzG1Wo1ff/0VCxcuxN27d/Hxxx8jJCQE/v7+YpdmsChcKkChUOCHH37A0qVLkZ6ejmHDhmHu3Llo1KiR2KXVaIwxPHnypNSBW3ONRskL/V4VOnXr1tVa6CgUilLXCmlqefjwoXCBaoMGDV4KES8vL7pAVWRqtRp79+7FwoULcf/+ffTp0wchISFo0aKF2KUZHAqXclAqldi2bRtCQ0ORmpqKIUOGYN68efDw8BC7NPIGPM/j0aNHL4VOXFycsLKunZ3dSwd7Hx8fODs7vzZ0lEqlsLROydd98OCBsLSOq6vrS6/p7e0Na2trnf39pPxUKhX27NmDRYsWISEhAZ9++ikWLFjw0koX5PUoXMqgqKgI27dvxzfffIOUlBQMHDgQ8+bNQ5MmTcQujVSCWq0W1sUqGQ537twRFld0dHSEt7c3XFxcYGlpCZVKhaysLNy/fx/x8fHCzavq1q37yhChRUENW3FxMcLCwrB48WI8evQI/fr1w4IFC+Dl5SV2aXqPwuUNiouL8dNPP2HJkiV48uQJBgwYgHnz5qFZs2Zil0aqgEqlwoMHDxAdHY0zZ84gKioK8fHxyMjIwIvNxNzcHG5ubvD19UWHDh3Qpk0b+Pj4wMHBQaTqSVUqKirCzp07Sx0L5s+fj6ZNm4pdmt6icHkFzbeVJUuWIDExEf369cP8+fPpLo/VRFnPWF4cl2nSpEmpa0w0z33bGYyPjw9sbW3F/JOJlrzYizFo0CDMmzcPnp6eYpemdyhcStD0sy5evBgPHjygflYDV9GxFl9fX9SuXbvMA/wlx15KbovGXqqvkuOvz549E8Zf3d3dxS5Nb1C44P9niCxatAjx8fHo06cPFixY8NKKwUQ/aWaJvTg1OS4uDgUFBQAAa2vrV4aINmeJvaiwsLDUbQ1o1lj18+LM0eHDh2Pu3Llo2LCh2KWJrkaHi2Zu+6JFi3Dnzh2a267nXnd9S2xsLPLz8wG8fH2L5qJKMa9vedHbrnfhOA6NGjWi610MyIvXvH3++eeYM2cO6tevL3ZpoqmR4aK5KnfhwoWIjY3FBx98gJCQELRu3Vrs0ghqzpX5L8rPz3/ptga3b98WrtSXSCTw8PCgK/X1WEFBAbZs2SKs1jFq1CjMnj27Rt52uUaFi2Y9oZCQENy+fRvvv/8+QkJC0LZtW7FLq7HS09NfGSJvW1OsYcOGNWYhyxfXGNOEztvWGGvSpAmtMC0SzTqDK1asgEwmw5gxYzBr1iy4uLiIXZrO1IhwoZVQxZebm4ubN2+WezVkd3d3g1/Jt6q8uDqyJnRKro7cpEmTUkvs+/r60vRZHdL2CumGpEaEy19//YWPPvqI7uEgoq1bt2LMmDF0HxcdeNN9XerVqyd0sxHdyc3Nxfr167F69Wo0btwYV69eFbukKlcjwkXzJ+rLgG5NVPJjRvtB9+j91w816VhUI8KFEEKIbulFZzZjDPHx8cLtWA2VRCKBr6+vQV6jQPtAfLQPxEf7QHv04syF53lMnDgRrq6usLKyEruccuF5HhkZGXBycsLZs2cxf/58g1yeuyL7QK1Wg+d5SKVSvTnNP3PmTLXYB8XFxTAxMTHIA3R12QdFRUWwsLCAmZmZ2GWVmz7sA704cwGeTzkdNWoUnJ2dxS6lzBhjOHPmDMaMGYNx48bB19f3pQUODUlZ9gFjDA8fPsTevXtx7tw5yGQyNGjQAD179kTPnj1hZ2cnWtAwxiCTyQx+H/Tq1QuzZs2Cs7Mz1q5dK+p7Wl7VZR98/PHHmDp1Ktzc3LB8+XI4OTnRPignvQkXQ5SVlYWZM2fi3r17OHLkCAICAsQuqUoVFRUhLCwMixYtwpMnT4R/P3v2LMLDw+Hv74/FixejW7duNeYalKpw/vx5nD17FoWFhVAqldi8eTPs7e0N5uBWHURGRuLixYuIiIjAgwcP8N1338HLy4v2QTkY5qXMekCtVmPNmjW4dOkSnJycsGTJkmq9NIdCocCSJUvwxRdfICkpCV5eXpg/fz42b96Mzz//HA4ODrh69Sr69++PtWvXQqlUil2ywerduzeWL18OCwsL/Pe//8WECROQk5Mj+jfRmqRPnz7YuHEjatWqhbNnz6JPnz6IjIwUFiIlb0fhUgGMMURGRmLjxo2QSCT4+uuv0apVK7HLqjLFxcVYtmwZli9fDrVajREjRuDkyZNYsGABxo4dix9++AEnTpxAr169IJPJMGfOHCxZskRYeZiUj1QqxdixY7Fy5UohYKZMmYL8/HwKGB0xNjbG0KFD8euvv6JJkyaIj4/HgAEDsG/fPuH2CuTNKFzKiTGGtLQ0zJw5E3l5eejWrRvGjRtXbU+XeZ7H9u3bsWLFCvA8j8mTJ2P9+vXC7X85jhNmpuzZswcTJkwAz/NYsWIFli1bhuLiYrH/BIMklUoxatQoLFu2DGZmZvj5558xffp0yOVyChgdkUgk6NKlC37//Xd07NgRaWlpGD16NLZs2UKf6zKgcCknlUqFpUuX4tq1a6hTpw6WLl1qcDPcykpzhjZnzhwolUoMHToUCxcuhIWFxUthynEcbGxssHz5ckyaNEkImK1bt9I3vQrSnMEsXLgQxsbG+PHHHxESEkJdjjrEcRyaNm2KX375Bf/6178gk8kwffp0LFu2jM7M34LCpRwYYzhy5Ai2bdsGIyMjzJkzBy1atKiWZy2MMTx69AhTpkxBVlYWunbtKowDvA7HcTA3N8eiRYswbNgwKJVKzJ07F0eOHKFv2xUklUoxefJkzJgxAxKJBOvXr8eqVavom7MOcRyHunXrYvv27Rg+fDiKi4uxZMkSzJs3DwUFBfTZfg0KlzLS3JBq1qxZkMvl+PjjjzFixAiDXd79bQoLCzFr1izExMSgQYMGWL9+PRwdHcsUpBYWFlixYgWCg4ORk5ODyZMnIyYmhhphBZmYmGDmzJmYOHEieJ5HaGgovvvuOzoj1CGO42BnZ4f169djypQpAIC1a9fi66+/prGw16ieR8YqUFRUhJCQEMTFxaFhw4YIDQ2ttrPDeJ7Htm3bsH//fpibm2P58uXw9vYu8xkax3Gwt7fHxo0b4e3tjcTEREyePBlZWVnUCCvI1NQUixYtwpAhQ6BUKjFnzhzs2bOHZi/pEMdxsLCwwJIlSzBr1ixIpVJs3boVkyZNotl8r0DhUgaMMfzyyy/Yu3cvTE1NsXjxYnh6elbb7rALFy5gyZIlUKvVGDduHPr06VPuv1VzN8UNGzbAwcEBp06dwpIlS6BSqaqo8uqN4zhYWlpi9erV6Nu3LwoKCjB16lT8/vvvFDA6xHEczMzMMHv2bCxevBimpqbYvXs3xo0bh8zMTAqYEihc3oIxhrt372LBggVQKpUYOHAgPv3002oZLACQlpaGadOmISMjAx07dsTs2bMrfD8VjuMQGBiIefPmQSqV4ocffsCBAweoAVaQpmtm06ZN6N69O7KzszFhwgQcO3aM3lMdMzExwZdffimMQ+7btw+jRo1Ceno67Yt/ULi8RWFhIebMmYNHjx7By8sLCxYsgImJidhlVZnTp08jJiYGzs7O+Pbbb+Hg4FCpIJVIJBgzZgw+/fRTFBYWYvbs2bh37x41wAriOA5OTk7YunWrMD121KhR+N///kfvqY5JpVKMGzcOa9asgbW1NQ4ePIjPP/8cqamptC9A4fJGPM9jx44dOHjwICwtLbFs2TK4ublV27MW4PnV4bt27cKaNWvQsmVLrfytZmZmCA0NhZeXFx4+fIiZM2dCLpdrodqaieM4uLq6YufOnWjTpg2ePn2KESNG4PTp03RQ0zGpVIoRI0Zgw4YNsLW1xeHDhzFs2DAkJyfX+H1B4fIajDHcunUL33zzDdRqNUaPHo2ePXtW62ABnl+Z3KdPH/Tv319rM+E4jkP9+vWxfPlyWFpa4tChQ9i6dSuNFVSCZkxr165d8Pf3R1JSEoYNG4YzZ87U+IOarhkZGWHQoEHCGnAnTpzA0KFD8eTJkxq9LyhcXkMmk2HWrFlITU2Fv78/ZsyYUWPu5a658l7br/n+++9j3LhxUKvVWLZsGa5du1ajG19lcRwHT09P7N69G++88w4eP36MoUOHUsCIQCKRoF+/fvjhhx/g6OiIyMhIDB48GImJiTV2X1C4vALP89i6dSuOHz8Oa2trLFu2DLVr1xa7LIMnlUoxffp0tGnTBmlpaZg+fTpyc3PFLsugcRyHZs2aYc+ePfDz88Pjx48xZMgQnDp1qsYe1MQikUjQp08f/Pjjj6hduzbOnj2LgQMH4sGDBzVyX1C4vIAxhhs3bghraU2YMAFBQUHVvjtMV2rVqoUVK1bA3t4ep0+fxvr16+liwErSBMzPP/+Md955B0+ePMGQIUNokF8EEokEH374IXbs2IE6derg0qVLGDhwYI2cxELh8gLNqr5paWlo3bo1pk2bRvcm0SKO49ChQwd8+eWXAJ5f5Xzu3Lka1/C0TRMwe/fuRatWrZCcnIyhQ4fi77//pvdWxyQSCd5//3389NNPcHFxwdWrVzFgwADExsbWqH1B4VKCpjvsxIkTsLGxwdKlS1GrVi2xy6p2jIyMMHnyZAQFBSEnJwfTp09HRkZGjWp4VYHjODRp0gQ///yzMIts+PDhOHToEE2e0DGO49C9e3eEhYWhfv36uHHjBj777DPcvHmzxnzOKVz+wRjDzZs3sXLlSvA8j/HjxyMwMJC6w6qIjY2NMJZ15coVrFixgrrHtIDjODRu3Bh79uxB+/bt8ezZM4wcORJ//vknBYyOcRyHrl27Yvfu3WjUqBFiYmLw2Wef4cqVKzUiYChc/lFQUIA5c+bg2bNnaNWqFXWHVTGO4+Dv74+ZM2fCyMgI33//PY4fP14jGl1V4zgO7u7u2LNnDzp37oz09HSMHj0aBw4coIDRMY7j0LFjR+zduxdNmjTBvXv3MHDgQJw/f77af9YpXPC8O+zHH38UZoeFhoZSd5gOSCQSjB49WriD5cyZM5GSklLtG50ucByHBg0aYPfu3QgKCkJmZibGjh2L8PBwOkPUMY7jEBAQgPDwcPj4+CAhIQGDBg1CZGRktf6s1/hwYYzh9u3bWL58OXiex9ixY9G1a1fqDtMRCwsLhIaGws3NDbdu3cLixYvpXiVaormSf9euXejRoweys7MxceJE7NmzhwJGxziOg5+fH8LDw4VrkoYMGVKtJ1zU+HCRy+WYM2eOcLHk119/XWMultQHHMeVWrNt165d+O2336ptg9M1juPg4uKCn376Cb169UJeXh4mTZqEn376iQJGxziOg7e3N8LDwxEQEICUlBQMHz4cf/31V7XsrqzR4cIYw44dO3Ds2DFYWVkhNDQUTk5OYpdV43AcJ6w2rVAoMHfuXCQkJFDAaAnHcXB2dsb27dvx8ccfQyaTYerUqfjhhx/oFgg6pllV4eeff0aHDh3w7NkzjBo1qlpOuKix4aLpDlu6dKmwdti7775L3WEiKXmfnAcPHmDu3Ll0j3ItKrmasuZ+MP/5z3+wadMm6obUMc2Ei927d6NLly7ChIvqdm+eGhsumu6wp0+fws/PD9OnT6fZYSLiOK7UHT4PHDiAsLAwOnvRIo7j4OjoiO+//x4DBgyAQqHA7NmzsWbNGhQVFYldXo2imXARFhYmTLgYM2ZMtZrRVyPDRbOU/tGjR4XusNq1a9NZi8g4jsPHH3+MESNGQKVSYdGiRbh9+zYFjBaVvAW15pbJCxYswLJly6BUKsUur0YpeeuEd999F1lZWRg3bly1OYOpceGiWUo/NDRU6A7r1q0bBYuekEqlmDNnDvz8/JCSkoJZs2ahoKBA7LKqFY7jYGtri3Xr1mHUqFFQqVQIDQ3F4sWLUVhYKHZ5NQrHcahXr16pgBk7diz++OMPgw+YGhcuJZfSb9myJWbMmEHdYXqE4zjUqVMHy5Ytg7W1NY4dO4bvv//e4BuavuE4DtbW1li9ejUmTpwInuexYsUKzJ07FwUFBXS2qEMlZ/S9++67wjVJhj6LrEaFC8/z2LJlC44fPy6sHebk5ERnLXqG4zi8++67GD9+vHDQo3u/aB/HcbC0tMTSpUsxbdo0cByH9evXY/r06ZDJZPR+65DmDGbHjh3o0qULMjIyMGbMGBw9etRg90ONCRfGGC5evCisHTZx4kRaSl+PSaVSfP3112jTpg3S09Mxa9Ys5Ofni11WtWRubo6QkBDMmjULUqkU33//PaZMmYLc3FyDPbAZopJjMB07dkRaWhpGjRqFEydOGOR+qBHhwhhDZmYm/vOf/yAzMxMdO3bEV199Rd1hes7R0RHLli2DnZ0dIiIisHnzZoPuJtBnZmZmmD17NkJCQmBqaopdu3Zh/PjxyMrKMsgDm6HS3BJ8165daNeuHVJTUzFixAiDXCpG78IlIyMDz5490+obqVKpsHTpUly8eBFOTk5YuXIl7O3ttfb6pGpwHIdOnTph0qRJAIBvv/0WV69eNbhGZihMTEwwbdo0YTr4r7/+itGjRyM9PZ3ecx3STMvftWsXWrduLVzJf/bsWYPaD3oVLrGxsejTpw9GjhyptdvfMsbw22+/4fvvv4eRkRHmzJmDgIAA6g4zEEZGRpgyZQratm2LjIwMzJ49G3l5eWKXVW0ZGxtj4sSJWL16NaysrPDHH39gxIgRSE1NNagDm6HjOA4eHh4ICwsT7i46bNgwXLp0yWD2g16FS0pKCm7fvo2jR49i8eLFlb6wizGGmJgYzJgxA3K5HJ9++ilGjRoFiUSv/mzyFvb29li6dCns7OwQGRmJLVu2UPdYFZJKpRg5ciQ2btwIW1tbHDlyBEOHDkVycrLBHNiqA83N38LCwuDr64vExEQMHToUUVFRBrEf9OooGxQUhLlz58LIyAhbtmxBeHh4hd9ExhiysrIwZcoUPHr0CD4+Pli2bBnMzc21XDWpapp7YnzxxRdgjOHbb781mAZmqIyMjDBo0CBs2bIFDg4OOHnyJIYOHYonT57Q+65DmsUuw8LC0KxZM9y/fx9Dhw5FdHS03u8HvQoXIyMjTJgwodTSFBWdgqpUKjFv3jxERkbCwcEB69atg5ubG3WHGShN91i7du2E7jGaPVa1JBIJ/v3vf+OHH35ArVq1EBkZicGDByMxMVHvD2zVCcdxaNGiBXbv3o3GjRvjzp07GDx4MGJiYvR6P+hVuADPZ62EhoaiZcuWSElJwZQpU8p9f3W1Wo0tW7Zg+/btMDY2RkhICN2jpRpwcHBAaGgobG1tERkZie+++466x6qYRCJB7969sW3bNjg5OeHs2bMYPHgwrVqtY5o7t4aFhQm3TB4yZAju3bunt/tB78KF4zjUrVsX69atg5OTEy5cuIB58+aVed0jnudx4MABhISEQKVSYcyYMRg9ejSNs1QDmtljX3zxBXiex+rVq+niSh2QSCT48MMPsX37djg7O+PChQsYOHAg7t+/T++9DnEchzZt2mDnzp1wc3PDzZs3MXToUDx8+FAv94NeHnE5jkP79u2xaNEimJiYYMeOHdiyZctbb27EGMOJEycwadIk5Ofn48MPPxReg1QPRkZG+PLLL9G2bVukp6dT95iOSCQS9OrVCz/99BNcXFxw5coVDBgwAHFxcXp5YKuuNOOPJffD8OHD9XKyhV6GC/D8wzx8+HCMGTMGarUaISEh+PXXX1/bDcLzPI4fP47PP/8caWlp6NixIzZt2gRbW1vqDqtmSnaPRURE0OwxHeE4Dj169MDOnTvh6uqKa9euoV+/fjS5Qsc4jkPXrl1LdVXq4/VIehsuwPOLuhYuXIgPP/wQ+fn5mDRpEsLDw0vdPY8xhqKiIuzevRtDhgxBSkoKWrdujR07dsDFxYWCpRriOA6dO3cu1T129epVscuqETTrvu3Zswfu7u6IjY1Fv379cOLECQp4HdIE/ebNm2FnZ4e///4bkyZN0qtrwPQ6XDRLg2/ZsgXdunUTlqOePXs27t69i/T0dFy8eBEjR47E2LFjkZ6ejs6dO+Pnn3+Gh4cHBUs1ZmRkhKlTp6J9+/bIyMjAjBkzkJOTI3ZZNYJm7OvXX39F8+bNkZiYiIEDB2Lnzp1022Qdkkgk6NOnD7799ltYWFhg//79mDVrlt7cwVUqdgFvo1mCfefOnZgyZQp+++03rFq1Ctu2bYOlpSWysrJQWFgIU1NTjBgxAt988w2cnZ0pWGoAe3t7LF++HP/6179w5swZrFu3jtaL0xHN7KV9+/ZhzJgxOH36NCZOnIh79+7B2NhY7PJqDIlEgsGDByMrKwtz5szBtm3b9ObGh3p95qKhCZgdO3Zg3bp18PHxgUKhwLNnz2BiYoKgoCD8/PPP2Lx5MwVLDaKZ+DFt2jRIJBKcO3cOMplM7LJqDI7j0LhxY/zyyy8YOnQo1Go14uPjqXtMx6RSKb744gtMnToVAHD58mW9uOmb3py5MMaQnZ391m89/fv3R8+ePZGQkICCggI4OzvDzc0NZmZmKCgoEPWuhYZ+m9iy7gN9M3DgQEilUvTu3Ru7d+8Wu5xKMcR9IJVKsWTJErRo0QJBQUE4cOCA2CVViiHuAwAYP348bGxs0Lt3b/z8889il6Mf4aJZZnrDhg0G3a1RWFgIW1tbscuoEEPfB4wxbNy4kfaBiBhjePjwIe0DETHGsHnzZr3YBxzTg7lrjDG9mkJXGRzHGWS3HO0D8dE+EB/tAy1uXx/ChRBCSPViEAP6lVWdvo0YKtoH4qN9IL6atA/0YsxFFzQ7VexTxZpG05BKNih6/8VRcl9QO9CtF9tBTXjva8SZCwAcPHgQLVu2hEQiQffu3XH+/HmxS6rW1Go19uzZg2bNmkEikaBv376Ijo6uEY1KXymVSqxfvx4uLi4wMTHB6NGjkZiYKHZZ1Vp+fj6++eYbODo6wtzcHF999RXS0tLELks3WA2iVqvZ/v37ma+vLwPAgoOD2cWLF8Uuq1pRq9Vs7969rFmzZgwA+/DDD9nVq1fFLouUUFBQwFatWsWcnJyYsbExGzt2LHv06JHYZVUr+fn5bOnSpczBwYGZmJiwSZMmseTkZLHL0qkaFS4aarWa/frrr8zb25sBYL169WJXrlwRuyyDpnlPfXx8GADWs2dPdunSJbHLIm8gk8nY8uXLWa1atZiJiQmbMGECe/LkidhlGTSZTMZWrFjBatWqxYyNjWv0e1ojw0VDpVKV+pb90UcfsaioKLHLMiias8HmzZsLZ4MXLlwQuyxSDiW/ZZuamtbIb9mVVVBQwFavXs1q167NpFIpnQ2yGh4uGiqViu3evZs1adKEAWC9e/dmN27cELssvcbzPPv999/ZO++8wwCwbt26sXPnzoldFqmE3NxctmTJEmZnZ8fMzMzYl19+yZ4+fSp2WXpNLpeztWvXsjp16jAjIyM2atQo9vDhQ7HL0gsULiUUFxeznTt3Mg8PDwaAffLJJyw6OlrssvQKz/Ps4MGDrFWrVgwACwoKYqdPnxa7LKJFOTk5bOHChczW1paZm5uzr776ij179kzssvRKYWEh27BhA6tbty4zMjJiI0aMYA8ePBC7LL1C4fIKxcXFbPv27axRo0YMAOvXrx+LiYkRuyxR8TzPDh8+zAICAhgA1rlzZ/a///1P7LJIFcrOzmbz589n1tbWzMLCgk2fPp2lp6eLXZaoFAoF27x5M3N1dWUSiYQNHTqUxcfHi12WXqJweYOioiK2detW1qBBA8ZxHBswYACLi4sTuyyd4nmeHTt2jLVr144BYB06dGAnTpxgPM+LXRrRkczMTDZnzhxmZWXFrKys2KxZs1hGRobYZemUUqlk33//PXNzc2Mcx7FBgwaxO3fuiF2WXqNwKQOlUsm+++475ubmxiQSCRs0aBC7e/eu2GVVKZ7n2YkTJ1iHDh0YANa2bVt27NgxCpUaLD09nc2cOZNZWloya2trNnfuXJaVlSV2WVWqqKiIbdu2TfiC+dlnn7HY2FixyzIIFC7loFAo2KZNm5iLi0u1PiWOiIhgXbp0YQBY69at2eHDhylUiCAtLY395z//Yebm5szGxoYtWLCAZWdni12WVhUXF7MdO3Ywd3d3BoD9+9//Zrdv3xa7LINC4VIBhYWFbP369cIMkeoymHf69GkWFBTEADB/f3928OBBChXyWqmpqWzatGnMzMyM2dnZsUWLFrHc3Fyxy6qU4uJitmvXLta4cWMGgPXt25fdvHlT7LIMEoVLJcjlcrZmzRrm7OzMpFIpGzVqFEtMTBS7rHI7d+4c69atGwPA/Pz82O+//06hQsosJSWFTZkyhZmamjJ7e3v2zTffsLy8PLHLKheVSsX27NkjXI7wr3/9i12/fl3ssgwahYsWvGo5jcePH4td1ltdvHiRBQcHMwDM19eX7d+/n6nVarHLIgYqKSmJffHFF8zExIQ5OjqyZcuWsfz8fLHLeiO1Ws3Cw8OZl5cXLVekZRQuWqRZTsPR0ZGZmJiwiRMnsqSkJLHLesmVK1dYr169GADm7e3Nfv31VwoVojWPHz9m48aNY8bGxqxWrVps5cqVrKCgQOyySlGr1ey///0vLVdUhShcqkBeXh4LDQ1l9vb2zNTUlE2ePJmlpKSIXRaLiopiH330EQPAmjZtyvbu3ctUKpXYZZFqKjExkY0ePZpJpVJWu3Zt9u233zK5XC5qTTzPswMHDrAWLVowAKxHjx7s/PnzotZUXVG4VKHc3Fy2ePFiYTmNqVOnstTUVJ3XcePGDda7d28GgHl6erLdu3dTqBCdSUhIYJ9//jkzMjJiderUYevWrWOFhYU6rYHnefbHH38IyxW999577OzZszqtoaahcNGB7OxsFhISIiyn8fXXX5drOQ2e51l6ejp7+PAhS09PL/Nge3R0NPvkk08YAObh4cF27tzJiouLK/pnEFIp9+/fZ8OHD2cSiYS5uLiwjRs3MoVCUebnV6Qd8DzP/vrrL2G5oq5du7JTp05V5s8gZUThokNZWVls3rx5wnIaM2bMeONyGtnZ2Wzt2rXCWmeaHw8PD7Z27drXXlsQExPD+vXrxwCwRo0ase3bt1OoEL1x7949NmTIECaRSJirqyvbsmULUyqVr/39irQDnufZkSNHWJs2bWi5IpFQuIggMzOTzZ49W1hOY/bs2SwzM7PU7xw9epRZWloyjuMYx3GlGpXm3ywtLdnRo0eF58TFxbEBAwYwjuNYgwYN2NatW1lRUZGu/zxCyuTOnTts4MCBjOM4Vr9+ffbDDz+89HktbzvgeZ79/fffpZYrOn78OE2tFwGFi4jS09PZjBkzmIWFBbO2tmbz5s1jWVlZ7OjRo8zIyIhJJJJSjenFH4lEwoyMjNiPP/7IBg8ezCQSCXNzc2PffffdG78JEqJPYmJiWP/+/RnHcaxRo0bsxx9/ZEVFReVqBxKJhC1btox17NhRWK7o6NGjFCoi4hhjTDs3TCYVlZaWhpUrV2LTpk0wNjZGYWEhVCoVyrNr6tati7lz52LkyJEwNTWtwmoJqRq3b9/GwoULsW/fPjRq1AjJyckoLi4uVzvw9/fHkiVL0LNnT3AcV4XVkrehcNEjqamp+Oyzz3Dq1KlyP3fVqlX46quvqqAqQnTr5s2bGDp0KKKjo8v93LVr12LKlClVUBUpLwoXPcIYg6enJx48eFCu53EcB3d3d8THx9O3NWLwqB1UDxQueiQjIwNOTk6Ver6jo6MWKyJE96gdVA8SsQsg/08mk1Xq+fn5+VqqhBDxUDuoHihc9IiVlVWlnm9tba2lSggRD7WD6oHCRY84OjrCw8Oj3P3FHMfBw8MDDg4OVVQZIbpD7aB6oHDRIxzHYdKkSRV67uTJk2kQk1QL1A6qBxrQ1zM5OTlwdXVFYWEheJ5/6+9LJBKYm5sjKSkJdnZ2VV8gITpA7cDw0ZmLnrGzs8P+/fvBcRwkkjfvHolEAo7jcODAAWpQpFqhdmD4KFz0UHBwMA4dOgRzc3NwHPfSab7m38zNzXH48GH06NFDpEoJqTrUDgwbhYueCg4ORlJSEtauXQt3d/dSj7m7u2Pt2rVITk6mBkWqNWoHhovGXAwAYwxZWVnIz8+HtbU1HBwcaNCS1DjUDgwLhQshhBCto24xQgghWkfhQgghROsoXAghhGgdhQshhBCto3AhhBCidRQuhBBCtI7ChRBCiNZRuBBCCNE6ChdCCCFaR+FCCCFE6yhcCCGEaB2FCyGEEK2jcCGEEKJ1FC6EEEK07v8Ae2oa6rMZtL0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 500x400 with 10 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "model.plot(beta=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "64d2573b",
   "metadata": {},
   "source": [
    "Fix the first layer activation to be linear function, and the second layer to be sine functions (caveat: this is quite sensitive to hypreparams)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "e2e78752",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "r2 is 0.8422812819480896\n",
      "r2 is not very high, please double check if you are choosing the correct symbolic function.\n",
      "r2 is 0.8454716801643372\n",
      "r2 is not very high, please double check if you are choosing the correct symbolic function.\n",
      "Best value at boundary.\n",
      "r2 is 0.9996306300163269\n",
      "r2 is 0.9994996190071106\n",
      "r2 is 0.998060405254364\n",
      "r2 is 0.9991359710693359\n"
     ]
    }
   ],
   "source": [
    "for i in range(2):\n",
    "    for j in range(2):\n",
    "        model.fix_symbolic(0,i,j,'x')\n",
    "        \n",
    "for i in range(2):\n",
    "    model.fix_symbolic(1,i,0,'sin')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3fae3f32",
   "metadata": {},
   "source": [
    "After setting all to be symbolic, we further train the model (affine parameters are still trainable). The model can now reach machine precision!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "308b72af",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "pde loss: 3.55e-11 | bc loss: 2.39e-13 | l2: 1.88e-13 : 100%|███████| 20/20 [00:11<00:00,  1.78it/s]\n"
     ]
    }
   ],
   "source": [
    "train()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35985ae9",
   "metadata": {},
   "source": [
    "Print out the symbolic formula"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "f0ec310e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - 0.5 \\sin{\\left(3.141592 x_{1} + 3.141593 x_{2} - 4.712389 \\right)} + 0.5 \\sin{\\left(3.141593 x_{1} - 3.141592 x_{2} + 1.570797 \\right)}$"
      ],
      "text/plain": [
       "-0.5*sin(3.141592*x_1 + 3.141593*x_2 - 4.712389) + 0.5*sin(3.141593*x_1 - 3.141592*x_2 + 1.570797)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "formula = model.symbolic_formula()[0][0]\n",
    "ex_round(formula,6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5c3e90ca",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
